source('../env.R')
Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
── Attaching core tidyverse packages ──────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.1 ✔ tibble 3.2.1
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
✔ purrr 1.0.2 ── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errorsRegistered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
Attaching package: ‘dbplyr’
The following objects are masked from ‘package:dplyr’:
ident, sql
Loading required package: ape
Attaching package: ‘ape’
The following object is masked from ‘package:dplyr’:
where
Loading required package: maps
Attaching package: ‘maps’
The following object is masked from ‘package:purrr’:
map
Please cite the eBird Status and Trends data using:
Fink, D., T. Auer, A. Johnston, M. Strimas-Mackey, S. Ligocki, O. Robinson,
W. Hochachka, L. Jaromczyk, C. Crowley, K. Dunham, A. Stillman, I. Davies,
A. Rodewald, V. Ruiz-Gutierrez, C. Wood. 2023.
eBird Status and Trends, Data Version: 2022; Released: 2023. Cornell Lab of
Ornithology, Ithaca, New York. https://doi.org/10.2173/ebirdst.2022
This version of the package provides access to the 2022 version of the eBird
Status and Trends Data Products. Access to the 2022 data will be provided
until November 2024 when it will be replaced by the 2023 data. At that
point, you will be required to update this R package and transition to using
the new data.
terra 1.8.5
Attaching package: ‘terra’
The following object is masked from ‘package:phytools’:
rescale
The following objects are masked from ‘package:ape’:
rotate, trans, zoom
The following object is masked from ‘package:tidyr’:
extract
Loading required package: permute
Loading required package: lattice
This is vegan 2.6-8
Attaching package: ‘vegan’
The following object is masked from ‘package:phytools’:
scores
Loading required package: mvtnorm
Loading required package: survival
Loading required package: TH.data
Loading required package: MASS
Attaching package: ‘MASS’
The following object is masked from ‘package:terra’:
area
The following object is masked from ‘package:dplyr’:
select
Attaching package: ‘TH.data’
The following object is masked from ‘package:MASS’:
geyser
Registered S3 method overwritten by 'broom':
method from
nobs.multinom MuMIn
Attaching package: ‘ggpubr’
The following object is masked from ‘package:terra’:
rotate
The following object is masked from ‘package:ape’:
rotate
Loading required package: nlme
Attaching package: ‘nlme’
The following object is masked from ‘package:dplyr’:
collapse
Registered S3 method overwritten by 'GGally':
method from
+.gg ggplot2
Attaching package: ‘GGally’
The following object is masked from ‘package:terra’:
wrap
Attaching package: ‘foreach’
The following objects are masked from ‘package:purrr’:
accumulate, when
Attaching package: ‘scales’
The following object is masked from ‘package:terra’:
rescale
The following object is masked from ‘package:phytools’:
rescale
The following object is masked from ‘package:purrr’:
discard
The following object is masked from ‘package:readr’:
col_factor
Attaching package: ‘cowplot’
The following object is masked from ‘package:ggpubr’:
get_legend
The following object is masked from ‘package:lubridate’:
stamp
── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────────────── treedataverse 0.0.1 ──
✔ tidytree 0.4.6 ✔ ggtree 3.14.0
✔ treeio 1.30.0 ✔ ggtreeExtra 1.16.0
── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────── treedataverse_conflicts() ──
✖ foreach::accumulate() masks purrr::accumulate()
✖ ggtree::collapse() masks nlme::collapse(), dplyr::collapse()
✖ scales::discard() masks purrr::discard()
✖ treeio::drop.tip() masks tidytree::drop.tip(), ape::drop.tip()
✖ ggtree::expand() masks tidyr::expand()
✖ tidytree::filter() masks dplyr::filter(), stats::filter()
✖ ggtree::flip() masks terra::flip()
✖ .GlobalEnv::geom_map() masks ggplot2::geom_map()
✖ dbplyr::ident() masks dplyr::ident()
✖ ggtree::inset() masks terra::inset()
✖ tidytree::keep.tip() masks ape::keep.tip()
✖ dplyr::lag() masks stats::lag()
✖ maps::map() masks purrr::map()
✖ treeio::mask() masks terra::mask()
✖ treeio::read.newick() masks phytools::read.newick()
✖ ggtree::rotate() masks ggpubr::rotate(), terra::rotate(), ape::rotate()
✖ tidytree::select() masks MASS::select(), dplyr::select()
✖ dbplyr::sql() masks dplyr::sql()
✖ terra::trans() masks ape::trans()
✖ foreach::when() masks purrr::when()
✖ ape::where() masks dplyr::where()
✖ terra::zoom() masks ape::zoom()
Loading required package: usethis
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Attaching package: ‘devtools’
The following object is masked from ‘package:permute’:
check
Using GitHub PAT from the git credential store.
Skipping install of 'clootl' from a github remote, the SHA1 (30ac57ee) has not changed since last install.
Use `force = TRUE` to force installation
To cite package ‘clootl’ in publications use:
Miller E, McTavish E, Sanchez-Reyes L (2025). “clootl: Fetch and Explore the Cornell Lab of Ornithology Open Tree of Life
Avian Phylogeny.” <https://github.com/eliotmiller/clootl>.
McTavish E, Gerbracht J, Holder M, Iliff M, Lepage D, Rasmussen P, Redelings B, Sanchez-Reyes L, Miller E (2025). “A complete
and dynamic tree of birds.” _Proceedings of the National Academy of Sciences_.
To see these entries in BibTeX format, use 'format(<citation>, bibtex=TRUE)', or 'toBibtex(.)'.
The current version of the Aves tree is v1.4.
Please specify the tree and taxonomy version used when citing this R package.
When possible, cite all the original studies supporting your tree:
These citations are acessible using getCitations(your_tree)
Loading required package: spData
To access larger datasets in this package, install the spDataLarge package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`
Registered S3 method overwritten by 'spdep':
method from
plot.mst ape
Attaching package: ‘gridExtra’
The following object is masked from ‘package:dplyr’:
combine
Attaching package: ‘grid’
The following object is masked from ‘package:terra’:
depth
Spherical geometry (s2) switched off
Load community data
communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))
Rows: 2428 Columns: 7── Column specification ────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, ebird_species_name, seasonal, presence, origin
dbl (2): city_id, relative_abundance_proxy
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
communities
Load trait data
traits = read_csv(filename(TAXONOMY_OUTPUT_DIR, 'traits_ebird.csv'))
Rows: 332 Columns: 10── Column specification ────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): ebird_species_name, habitat, trophic_level, trophic_niche, primary_lifestyle
dbl (5): beak_width, hwi, mass, habitat_density, migration
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(traits)
species_urban_response = communities %>% group_by(ebird_species_name) %>% summarise(mean_relative_abundance_proxy = mean(relative_abundance_proxy)) %>% left_join(traits)
Joining with `by = join_by(ebird_species_name)`
species_urban_response$is_urban = species_urban_response$mean_relative_abundance_proxy > 0
species_urban_response$ebird_species_name = str_replace(species_urban_response$ebird_species_name,' ', '_')
species_urban_response
species_urban_response %>% filter(mean_relative_abundance_proxy > 50)
species_urban_response %>% filter(mean_relative_abundance_proxy > 40)
species_urban_response %>% filter(mean_relative_abundance_proxy > 30)
species_urban_response %>% group_by(is_urban) %>% summarise(
median_beak_width = median(beak_width),
iqr_beak_width = IQR(beak_width),
median_hwi = median(hwi),
iqr_hwi = IQR(hwi),
median_mass = median(mass),
iqr_mass = IQR(mass),
n = n()
)
ggplot(species_urban_response, aes(x = is_urban, y = beak_width)) + geom_boxplot()
wilcox.test(beak_width ~ is_urban, species_urban_response)
Wilcoxon rank sum test with continuity correction
data: beak_width by is_urban
W = 1525, p-value = 0.3893
alternative hypothesis: true location shift is not equal to 0
There is not a significant response for the gape width between urban species and non urban species.
ggplot(species_urban_response, aes(x = is_urban, y = hwi)) + geom_boxplot()
wilcox.test(hwi ~ is_urban, species_urban_response)
Wilcoxon rank sum test with continuity correction
data: hwi by is_urban
W = 1213, p-value = 0.01492
alternative hypothesis: true location shift is not equal to 0
There is a significant response between urban and non-urban species
ggplot(species_urban_response, aes(x = is_urban, y = mass)) + geom_boxplot()
wilcox.test(mass ~ is_urban, species_urban_response)
Wilcoxon rank sum test with continuity correction
data: mass by is_urban
W = 1825, p-value = 0.5168
alternative hypothesis: true location shift is not equal to 0
There is not a significant response between species.
tree = read.tree(filename(TAXONOMY_OUTPUT_DIR, 'phylogeny.tre'))
tree_pruned = ladderize(drop.tip(tree, setdiff(tree$tip.label, species_urban_response$ebird_species_name)))
ggtree(tree_pruned, layout="fan")
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
columbidae_response <- species_urban_response$mean_relative_abundance_proxy
names(columbidae_response) <- species_urban_response$ebird_species_name
Is the a phylogenetic signal for our relative abundance proxy
col_ut_phylo_signal <- phylosig(tree_pruned, columbidae_response, method="lambda", test=TRUE)
col_ut_phylo_signal
Phylogenetic signal lambda : 0.46832
logL(lambda) : -576.839
LR(lambda=0) : 7.45324
P-value (based on LR test) : 0.00633224
species_urban_response$ebird_species_name_clean = gsub("_", " ", species_urban_response$ebird_species_name)
normalize <- function(x) (x - min(x)) / (max(x) - min(x)) * 5 + 1 # Scale to range [1, 6]
species_urban_response$beak_width_norm <- normalize(species_urban_response$beak_width)
species_urban_response$hwi_norm <- normalize(species_urban_response$hwi)
species_urban_response$mass_norm <- normalize(species_urban_response$mass)
g = ggtree(tree_pruned, layout="fan") %<+% species_urban_response + abundance_proxy_scale
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
g +
geom_fruit(geom = geom_point, mapping = aes(size = beak_width_norm), fill = beak_width_colour, colour = 'black', shape = 21, stroke = 0.2, offset = 0, show.legend = FALSE) +
geom_fruit(geom = geom_point, mapping = aes(size = hwi_norm), fill = hwi_colour, colour = 'black', shape = 21, stroke = 0.2, offset = 0.13, show.legend = FALSE) +
geom_fruit(geom = geom_point, mapping = aes(size = mass_norm), fill = mass_colour, colour = 'black', shape = 21, stroke = 0.2, offset = 0.15, show.legend = FALSE) +
geom_tiplab(size=2, aes(label=ebird_species_name_clean, color = mean_relative_abundance_proxy), offset = 15) +
labs(color = "Mean relative abundance proxy") +
theme(legend.position = "bottom") + ggplot2::xlim(0, 80)
ggsave(filename(FIGURES_OUTPUT_DIR, 'phylogeny.jpg'), width = 2000, height = 2100, units = 'px')
Format plot for nature
colors = c('Beak Width' = beak_width_colour, 'HWI' = hwi_colour, 'Mass' = mass_colour)
g_fig_base = ggtree(tree_pruned, layout="fan") %<+% species_urban_response
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
g_fig = g_fig_base +
geom_fruit(geom = geom_point, mapping = aes(size = beak_width_norm, fill = 'Beak Width'), colour = 'black', shape = 21, stroke = 0.2, offset = 0, show.legend = c(size = FALSE, fill = TRUE)) +
geom_fruit(geom = geom_point, mapping = aes(size = hwi_norm, fill = 'HWI'), colour = 'black', shape = 21, stroke = 0.2, offset = 0.13, show.legend = c(size = FALSE, fill = TRUE)) +
geom_fruit(geom = geom_point, mapping = aes(size = mass_norm, fill = 'Mass'), colour = 'black', shape = 21, stroke = 0.2, offset = 0.15, show.legend = c(size = FALSE, fill = TRUE)) +
geom_tiplab(size=2, aes(label=ebird_species_name_clean), offset = 30) +
scale_fill_manual(name = "Relative trait size", values = colors) +
guides(fill = guide_legend(override.aes = list(size = 6))) + new_scale_fill()
heatmap_fig_ata = data.frame(species_urban_response[, "mean_relative_abundance_proxy"])
colnames(heatmap_fig_ata) = c('MRAP')
rownames(heatmap_fig_ata) = species_urban_response$ebird_species_name
gheatmap(g_fig,
heatmap_fig_ata,
offset = 14,
width = 0.2,
colnames_angle = 0,
colnames_offset_x = 0.5,
font.size = 3) +
scale_fill_continuous(low = cyan, high = orange, guide = guide_colorbar(barwidth = 11)) +
labs(fill = "MRAP\n(Mean relative abundance proxy)") +
theme(legend.position = "bottom", legend.direction = "horizontal", legend.box = "vertical", legend.title.position = "top", legend.justification = "centre") +
ggplot2::xlim(0, 85)
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
ggsave(filename(FIGURES_OUTPUT_DIR, 'figure3.jpg'), width = 180, height = 220, units = 'mm', dpi = 450)
cut_off_3 = 50
cut_off_2 = 35
cut_off_1 = 1
hwi_by_mass_urban_3 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_3) %>% slice(chull(hwi, mass))
hwi_by_mass_urban_2 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_2) %>% slice(chull(hwi, mass))
hwi_by_mass_urban_1 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_1) %>% slice(chull(hwi, mass))
beak_width_by_mass_urban_3 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_3) %>% slice(chull(beak_width, mass))
beak_width_by_mass_urban_2 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_2) %>% slice(chull(beak_width, mass))
beak_width_by_mass_urban_1 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_1) %>% slice(chull(beak_width, mass))
hwi_by_beak_width_urban_3 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_3) %>% slice(chull(hwi, beak_width))
hwi_by_beak_width_urban_2 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_2) %>% slice(chull(hwi, beak_width))
hwi_by_beak_width_urban_1 = species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_1) %>% slice(chull(hwi, beak_width))
species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_3) %>% nrow()
species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_2) %>% nrow()
species_urban_response %>% filter(mean_relative_abundance_proxy > cut_off_1) %>% nrow()
polygon_line_type = 'dashed'
polygon_linewidth = 0.4
polygon_alpha1 = 0.05
polygon_alpha2 = 0.1
polygon_alpha3 = 0.15
hwi_by_mass = ggplot(species_urban_response, aes(x = hwi, y = mass, colour = mean_relative_abundance_proxy)) +
geom_polygon(data = hwi_by_mass_urban_1, alpha = polygon_alpha1, color = trait_polygon_light, fill = trait_polygon_light, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_polygon(data = hwi_by_mass_urban_2, alpha = polygon_alpha2, color = trait_polygon_med, fill = trait_polygon_med, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_polygon(data = hwi_by_mass_urban_3, alpha = polygon_alpha3, color = trait_polygon_dark, fill = trait_polygon_dark, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_point() +
abundance_proxy_scale +
labs(colour = "Mean relative abundance proxy", y = 'Mass', x = '') +
theme(legend.position = 'bottom', legend.title.position = 'top', legend.key.width = unit(12, 'mm'))
beak_width_by_mass = ggplot(species_urban_response, aes(x = beak_width, y = mass, colour = mean_relative_abundance_proxy)) +
geom_polygon(data = beak_width_by_mass_urban_1, alpha = polygon_alpha1, color = trait_polygon_light, fill = trait_polygon_light, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_polygon(data = beak_width_by_mass_urban_2, alpha = polygon_alpha2, color = trait_polygon_med, fill = trait_polygon_med, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_polygon(data = beak_width_by_mass_urban_3, alpha = polygon_alpha3, color = trait_polygon_dark, fill = trait_polygon_dark, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_point() +
abundance_proxy_scale +
labs(y = '', x = 'Beak width')
hwi_by_beak_width = ggplot(species_urban_response, aes(x = hwi, y = beak_width, colour = mean_relative_abundance_proxy)) +
geom_polygon(data = hwi_by_beak_width_urban_1, alpha = polygon_alpha1, color = trait_polygon_light, fill = trait_polygon_light, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_polygon(data = hwi_by_beak_width_urban_2, alpha = polygon_alpha2, color = trait_polygon_med, fill = trait_polygon_med, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_polygon(data = hwi_by_beak_width_urban_3, alpha = polygon_alpha3, color = trait_polygon_dark, fill = trait_polygon_dark, linewidth = polygon_linewidth, linetype = polygon_line_type) +
geom_point() +
scale_fill_continuous(low = cyan, high = orange, guide = guide_colorbar(barwidth = 100)) +
labs(y = 'Beak width', x = 'HWI')
legend <- ggpubr::get_legend(
# create some space to the left of the legend
hwi_by_mass
)
plot_grid(
nrow = 2, ncol = 2,
hwi_by_mass + theme(legend.position="none"),
beak_width_by_mass + theme(legend.position="none"),
hwi_by_beak_width + theme(legend.position="none"),
legend
)
ggsave(filename(FIGURES_OUTPUT_DIR, 'traits.jpg'), width = 2000, height = 1800, units = 'px')